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ON THE DYNAMICS OF SOME GRID ADAPTION SCHEMES* 1 ') 
P.K.Sweby*”) and H.C. Yee<™> 


The dynamics of a one-parameter family of mesh eqnidistribution 
schemes coupled with finite difference discretisations of linear and non- 
linear convection- diffusion model equations is studied numerically. It is 
shown that, when time marched to steady state, the grid adaption not only 
influences the stability and convergence rate of the overall scheme, but can 
also introduce spurious dynamics to the numerical solution procedure. 

1. INTRODUCTION 

It has been shown recently by the authors and others (see e.g. [1-4] 
and references therein), that the dynamics of the numerical discretisations 
of non-linear differential equations (DEs) can differ significantly from that 
of the original DEs themselves. For example, the discretisations can pos- 
sess spurious numerical steady- state solutions and spurious high order 
period solutions (oscillatory behaviour) which may be stable or unstable 
(but still affecting the allowable initial data which will converge to the 
stable steady-state solutions) and can occur below or above the linearised 
stability limit on the time step for the numerical scheme. These studies 
are particularly important for computational fluid dynamics (CFD) since 
it is common practice to use a time- dependent approach to obtain steady 
state numerical solutions of complicated steady fluid flows. References [1- 
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4] have demonstrated, using simple test cases, how such computations may 
converge to incorrect solutions, be non- convergent or suffer from residual 
plateauing. 

Grid adaption is commonly used in CFD applications to economically 
resolve strong gradients and complex shock wave interactions. It is well 
known (see e.g. [5]) that some grid adaption techniques for time-dependent 
partial differential equations (PDEs) can experience instability, depending 
on the methods used in solving these grid adaption equations. 

In this paper we apply techniques used in [1-4] to study the dynamics 
of grid adaption using the equidistribution principle for time- dependent 
approach to steady-state solutions of problems involving steep gradients. 
In a parallel study [6] use is made of the AUTO computer bifurcation 
package [7] to obtain bifurcation diagrams for similar grid adaption meth- 
ods for the steady PDEs. However, the dependence on known solutions 
of the discretised PDEs and grid equations as starting values limits its 
usage. In the present work we utilise the power of the highly parallel Con- 
nection Machine CM-5 to undertake a purely numerical investigation into 
the dynamics possible in grid adaption schemes for the time-dependent 
PDEs. 

Section 2 discusses the equidistribution technique and weight func- 
tion used. Section 3 describes the test problems whilst Section 4 presents 
selected results and some concluding observations. 

2. GRID ADAPTION 

One common criterion used for grid adaption is the equidistribu- 
tion of a positive definite weight function u>(x,t), often taken to be some 
monitor of the numerical solution u(*,t) of the underlying PDE. A grid 
*o < *i(t) < ••• < xjv_i(t) < *jv, where Zf> and zjf are fixed, equidis- 
tributes w(z,t ) (at time t ) if 

/ * w(z,t)dz = f k+ ' w(z,t)dz = -^ / " w{z,t)dz, (1) 

for k = 1, . . . N. Here we choose a one-parameter family of weight func- 
tions 

w(z,t) = y/l - a + auj(*,f), a 6 [0,1], (2) 

where a = 1/2 corresponds to equidistribution in the arc-length and a = 0 
yields a uniform grid. 

If we now approximate w(z, t) to be constant in each interval (**_!, z h ) 



we have 


( Jl - a + au* ) (* fc - **_! ) = ( Jl- a + auj (x fc+ i - ®j,). 

\ V / h- 1/2 v ' k+1/2 

(3) 

Thus given a numerical solution of the PDE we can approximate the 
derivatives by 

I U * ~ lA\ 

U.U-1/2 - • ( 4 ) 

- * l »- 1 


Equation (3) is nonlinear in {«*} using (4). However, (3) is linear in { xj , } if 
{**} in (4) uses the existing grid. In this case we can solve the tridiagonal 
system (3) for a new set of { x to obtain an updated grid, the modified 
solution values on which may be obtained either by interpolation from the 
old grid or a convective update to take into account the grid movement. 
The majority of our computations employ this linearization of {**} in (4). 
(Note that [6] solves a slightly different form than (3).) 


Various strategies are possible using the regriding technique (3). For 
example, the grid adaption can be applied after every time step of solving 
the PDE or after a prescribed number of time step based on some criteria 
of the computed solutions. The {u*} can be updated via interpolation 
from the old values. If time-marching to the steady-state numerical so- 
lution is sought, crude approximations to updating solution values {u*} 
immediately following a regriding step may be made. For example inter- 
polation from the old grid can be used even after large grid movement or 
no adjustment of solution values {tx*} at all. In such cases we let the {tx*} 
be adjusted by the numerical scheme itself as it converges toward a steady 
state. 


3. TEST PROBLEMS 

The test problems used in this study are linear and non-linear forms 
of the convection-diffusion equation 

«t + /(«). = eu„ (5) 

with /(u) = u and f(u) = ^ti 2 respectively. We impose boundary condi- 
tions u(0,t) = 0 and u(l,t) = 1 for the linear case and u(0,t) = 1 and 
u(l,t) = —1 for the non-linear case which result in steady-state solutions 
of a boundary layer at x = 1 and a viscous shock at x = 1/2 respectively. 
In both cases the steepness of the feature is governed by the parameter e. 

We use a method of lines approach to solve the PDE with central 
spatial differencing for the diffusion term and either upwind or central 
differencing for the convective term. The resulting system of ordinary 



differential equations (ODEs) for du h /dt is then solved using a standard 
numerical method for ODEs. We have experimented with various explicit 
methods, such as forward Euler, second and fourth order Runge-Kutta, as 
well as the linearised implicit (backward) Euler method. As expected the 
stable time step required for the explicit methods was orders of magnitude 
lower than that for the implicit method and so we concentrated on this 
latter method of solution in this paper. Note also that unlike the stan- 
dard explicit Runge-Kutta methods (i.e. all but forward Euler), linearised 
implicit Euler cannot introduce spurious steady states. However, it can 
stabilise genuine unstable steady states of the system; see for exam ple 
[2]. Note that for this scheme and the linear problem, any non-linearity 
present is due solely to the grid adaption. 

The regriding strategy adopted was to regrid after every time step of 
the PDE method, either interpolating updated solution values Rom the 
old grid or performing no adjustment at all due to grid movement. This 
latter approach in effect presents the PDE method with new initial data 
to the problem at each step. 

We performed many of our calculations on a CM-5 connection ma- 
c hine at NASA Ames. Its parallel architecture enabled us to easily perform 
computations covering a range of selected parameters, usually e, but also 
the time step At and the monitor parameter a. We use nodal placement 
and the l 2 norm of the solution to illustrate our results. In all of the 
computation on the CM-5, we use the previous time step value for sc* in 
(4) to achieve a linear tridiagonal system for the updated grid in (3). 

For the majority of the results, 19 free (interior) nodes were used. 
Numerical experiments indicate that there is no major difference in quality 
of the results with either an odd or an even number of nodes. We divide 
a chosen parameter space (e.g., e) into 512 equal intervals with the rest of 
parameters (a, At , initial data) fixed. For each chosen parameter value, we 
iterate the discretised PDE and the grid function 4,000 steps (8,000 steps 
for the explicit methods) to allow the solution to settle to an asymptotic 
state. Then, we perform a series of time step /regriding stages, during 
which we investigate the dynamics by producing an overlaid plot of the 
norms at each step, resulting in a bifurcation type diagram. 

To examine the effect of grid adaption alone compared with the time- 
dependent approach to the steady states with regriding above, we also 
solve the regriding equations (3) and (4) iteratively using the analytic 
solutions of the PDE (5), starting from a uniform grid. We also used a 
quadrature based equidistribution applied to (1), namely the trapezium 



rule with a large number of subdivisions, to obtain the “exact” placements 
of the grid nodes for the analytical steady-state solution. 

Due to a page limitation, details of the formulas, procedures and 
extended results on the above will be reported in an extended version of 
this paper. The following summarized a small portion of the results using 
central difference for the convection and diffusion terms and the linearized 
implicit Euler method. 

4. RESULTS 

Although results for various values of a were computed, we focus on 
only a few of those obtained using a = 1/2, but will remark on results 
obtained with other values towards the end of this section. 

4.1 The Linear Problem. 


Figure 1 shows the nodal position for the steady-state solution ob- 
tained using our regriding technique (3) iteratively on the analytic form 
of the solution of the PDE (5). As expected the grid is concentrated 
around x = 1 for small e due to the very steep gradients of the solution 
there. These results agree with the “exact” nodal placements obtained 
by the fine quadrature equidistribution mentioned in the previous section. 
The figure shows a settled solution, which took between 100 iterations for 
the smallest value e and 10 iterations for the largest. It is interesting to 
note that a staircase effect is observed on the plots with iterations against 
logic e (figures not shown). 

Figure 2 shows the corresponding nodal positions when the time de- 
pendent PDE is solved to steady state using central spatial differences, 
with no adjustment to {u*} after regriding. Unless we state otherwise 
the nodal position plots show the final 5 nodal positions for each node, 
overlaid on the same plot. In this particular case it can be seen that for 
e > 10“ 2 6 the grid has settled, whilst for lower values of e this is not the 
case. Notice the jump in nodal positions as e is varied through 10“ 2,4 . 
This is a phenomenon observed in many of our calculations, although the 
particular value of e at which this would happen varied. Careful examina- 
tion near e = 10“ 3,5 suggests some slight structure in the pattern of nodal 
positions. This is revealed more clearly in Figure 3, a bifurcation diagram 
of the ( 2 norm of the solution, where a structure resembling a period seven 
bifurcation can be seen. For larger values of e the norm of the numerical 
solution agrees with that of the analytic steady state. 
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Figure 3: Solution norm for the linear problem, central differences, At - 1 

Figure 4 shows the bifurcation diagram for a selected range of At 
values. Note that the selected range of At is divided into 512 equal inter- 
vals. For each (fixed) At value, a fresh calculation of the discretised PDE 
and grid function was performed The correct solution is obtained for both 
large and small values of At, but additional spurious asymptotes occur for 
the range 0.1 < At < 10. In our other work [2] we have found that this 
is not un-typical of the linearised implicit Euler scheme used here. The 
nodal positions are given in Figure 5 and agree reasonably well, where 
they have converged, with other results using the same value of e. It is to 
be noted that converged nodal positions are independent of the time step 
At of the scheme. 

The previous results were produced without adjusting nodal values 
between regriding and the subsequent time step iteration (a procedure 
which can only be contemplated for situations where a steady state is 
sought because the transient behaviour of the solution is unimportant). 
Figures 6 and 7 show results for the same parameter values, but where 
interpolation using the old grid has been used to adjust the solution values 
at each regriding step. As can be seen this has modified the behaviour but 
not completely eliminated the problem. Notice in particular how there is 
a greater amount of node movement compared to the norm of the solution 
than in the previous case. Also note that the range of At effected has 
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Figure 4: Solution norm for the linear problem, central differences e = 

0.001 



Figure 5: Last 5 nodal positions for the linear problem, central differences 
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Figure 6: As Figure 4 but using interpolation 


shifted. 

4.2 The Non-linear Problem. 

Figure 8 shows the final 5 sets of nodal positions obtained by apply- 
ing the regriding technique iteratively to the analytic steady-state solu- 
tion. Here the nodal positions are incorrect for small e when compared 
with those obtained using quadrature based equidistribution. The correct 
positions can be visualised by extrapolating back from those for e = 10 -1 . 
Notice again the jumps in the nodal positions as e is varied. In the middle 
range of e a period two solution can be observed for the centre three nodes. 
This is one case where the effect is much more pronounced for an even 
number of free nodes. 

Much of the dynamics observed in the linear problem can be found in 
the non-linear case, although due to its increased severity the features are 
often more pronounced and convergence less readily attained. Detailed 
results will be reported in an extended version of this paper. To give a 
flavor of the dynamics for the nonlinear case, Figure 9, shows an overlaid 
plot of the values taken by the numerical solution itself for one set of 
parameter values in the range where convergence is not obtained. 
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Figure 8: Last 5 nodal positions for the exact non-linear problem using 
matrix iteration 
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Figure 9: Superposition of last 512 solution profiles for non-linear problem, 
centred differences, A t = 0.01, e = 0.01 

4.3 The effect of a. 

As mentioned previously we have also carried out calculations with 
other values of a in our weight function (2). A value of a = 0 yields a 
uniform gird, whilst a = 1 is in general unsuitable since strict positivity 
of the weight function is lost. This is particularly the case in our test 
problems for small e due to the solution gradients being very near zero. For 
a values other than 1, similar phenomena are encountered as for a = 1/2. 
However, the nodes are more tightly clustered around the steep gradients 
of the solution for larger a and are more spaced out for lower a. Also, for 
a near 1, convergence becomes harder to obtain. There is an increase in 
spurious dynamics of the numerical solutions and nodal values, especially 
for the non-linear problem. Figure 10 illustrates some of the behaviour 
for the linear problem. 

4.4 Concluding Remarks. 

We have presented a small selection of results from our preliminary 
investigation into the dynamical properties of grid adaption schemes. Fu- 
ture research will expand on the present work, and in particular concen- 
trate on understanding the complicated dynamical behaviour in the hope 




Figure 10: Dependence on a for the linear problem using interpolation, 
central differences, At = 1, c = 0-01 

of mapping out a reliable range of grid weight functions and grid adap- 
tion parameters of commonly used numerical schemes in CFD for flows 
containing strong gradients. 

REFERENCES 

1 YEE, H.C., SWEBY, P.K., GRIFFITHS, D.F. — Dynamical Ap- 
proach Study of Spurious Steady- State Numerical Solutions for Non- 
linear Differential Equations, Part I, The Dynamics of Time Dis- 
cretizations and its Implications for Algorithm Development in Com- 
putational Fluid Dynamics. J. Comput. Phys. , Vol. 97, pp 249-310, 
1991. 

2 YEE, H.C., SWEBY, P.K. — Dynamical Approach Study of Spu- 
rious Steady-State Numerical Solutions for Nonlinear Differential 
Equations, Part II, The Dynamics of Numerics of Systems of 2 x 
2 ODEs and its Connections to Finite Discretizations of PDEs. 
RNR Technical Report RNR-92-008 , March 1992. 

3 LAFON, A., YEE, H.C. — Dynamical Approach Study of Spurious 
Steady-State Numerical Solutions of Nonlinear Differential Equa- 
tions, Part III, The Effects of Nonlinear Source Terms and Bound- 


ary Conditions in Reaction- Convection Equations. NASA Technical 
Memorandum 103877 , July. 1991. 

4 LAFON, A., YEE, H.C. — Dynamical Approach Study of Spu- 
rious Steady- State Numerical Solutions of Nonlinear Differential 
Equations, Part IV, Stability vs. Numerical Treatment of Nonlinear 
Source Terms. ONERA/CERT Technical Memorandum , Feb. 1992. 

5 COYLE, J.M., FLAHERTY, J.E. & LUDWIG, R. — On the Stabil- 
ity of Mesh Equi distribution Strategies for Time-Dependent Partial 
Differential Equations. J. Comput. Phys. , Vol. 62, pp26-39, 1986. 

6 BUDD, C.J., STUART, A.M., KOOMULLIL, G.P. YEE, H.C. — 
Numerical Solution Behaviour of Model Convection— Diffusion BVP 
with Grid Adaptation. In preparation . 

7 DOEDEL, E. — AUTO: Software for Continuation and Bifurcation 
Problems in Ordinary Differential Equations, CIT Press, Pasadena, 


1986. 



